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ABSTRACT 

Intergalactic filaments form the foundation of the cosmic web that connect galaxies together, and 
provide an important reservoir of gas for galaxy growth and accretion. Here we present very high 
resolution two-dimensional simulations of the thermal and chemical evolution of such filaments, making 
use of a 32 species chemistry network that tracks the evolution of key molecules formed from hydrogen, 
oxygen, and carbon. We study the evolution of filaments over a wide range of parameters including 
the initial density, initial temperature, strength of the dissociating UV background, and metallicity. 
In low-redshift, Z rs Q.1Z@ filaments, the evolution is determined completely by the initial cooling 
time. If this is sufficiently short, the center of the filament always collapses to form dense, cold core 
containing a substantial fraction of molecules. In high-redshift, Z = 10~ 3 Zq filaments, the collapse 
proceeds much more slowly. This is due mostly to the lower initial temperatures, which leads to a 
much more modest increase in density before the atomic cooling limit is reached, making subsequent 
molecular cooling much less efficient. 

Finally, we study how the gravitational potential from a nearby dwarf galaxy affects the collapse of 
the filament and compare this to NGC 5253, a nearby starbusting dwarf galaxy thought to be fueled 
by the accretion of filament gas. In contrast to our fiducial case, a substantial density peak forms at 
the center of the potential. This peak evolves faster than the rest of the filament due to the increased 
rate at which chemical species form and cooling occur. We find that we achieve similar accretion rates 
as NGC 5253 but our two-dimensional simulations do not recover the formation of the giant molecular 
clouds that arc seen in radio observations. 



1. INTRODUCTION 

In the modern theory of cosmic evolution, structures 
formed from the bottom up through continuous accretion 
and mergers (e.g. White & Rees 1978; White & Frenk 
1991; Kauffmann et al 1993; Cole et al 2000; Bower et al 
2006). Originating as small perturbations in the initial 
density field, seed overdensities were slowly transformed 
by non-linear growth into the complex distribution of 
galaxies observed in large scale surveys (e.g. Doroshkc- 
vich et al 2004; Ratcliffc et al 1996; Pandey & Bharadwaj 
2006) and seen in TV-body simulations (e.g. Bond, Kof- 
man, & Pogosyan 1996). Between these growing struc- 
tures is a complex network of filaments, which connects 
them together and plays a key factor in their evolution. 

This cosmic web provides a reservoir of gas by which 
galaxies grow and accrete (White & Rees 1978), which 
may explain the nearly constant star-formation experi- 
enced by the Milky Way for example (Binney, Dehnen 
& Bertelli 2000). The state of the gas flowing along a 
filament is predominately determined by the mass of the 
dark matter halo (e.g. Birnboim & Dekel 2003; Dekel 
& Birnboim 2006; Keres et al 2005). Traditionally, two 
modes of accretion are discussed. "Hot mode" accretion 
occurs when the inflowing gas is shock heated as it is com- 
pressed by the hot gas of the halo and typically occurs 
in galaxies with masses above 10 12 M Q (Thoul & Wein- 
berg 1995; Furlanetto & Loeb 2004; Birnboim & Dekel 
2003; Dekel & Birnboim 2006). Below this threshold, 
the gas is not shock heated and accretes directly onto 
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the galaxy. At high redshifts, the cold mode dominates 
as most galaxies are not massive enough to shock heat 
the infalling gas. Conversely, the hot mode dominates at 
low redshift (z < 2) where most massive galaxies form 
(Kacprzak et al 2011; Stewart et al 2011). 

In addition, filamentary accretion may help explain the 
formation of the observed cold neutral hydrogen clouds 
that surround the Milky Way. These high- velocity clouds 
(HVCs) could form in the Milky Way halo from a vari- 
ety of mechanisms: cooling instabilities in the hot cir- 
cumgalactic medium (Mo & Miralda-Ecude 1996; Mailer 
& Bullock 2004; Joung et al 2012), production from a 
galactic fountain where hot gas is driven from the disk, 
cools and condenses, and rains back onto the disk (Breg- 
man 1980; Melioli et al 2008, 2009; Marinacii et al 2011), 
or the interaction between hot halo gas and an infalling 
satellite (e.g. Putman et al 2003). On the other hand, 
Keres et al (2009) suggest that through turbulent frag- 
mentation, the filament gas itself is able to create these 
cold clouds. 

Observing filaments has long been a challenging task. 
Direct detection attempts initially concentrated on X- 
ray surveys, motivated by the fact that most of the fil- 
amentary gas is expected to be very hot (e.g. Pierre, 
Bryan, & Gastaud 2000). Only recently have such stud- 
ies been successful (Werner et al 2008; Walker et al 2012; 
Williams et al 2010). Weak lensing, on the other hand, 
have has some success in studying filaments of galaxy 
cluster pairs (Gray et al 2002; Dietrich et al 2005). Still 
others have derived filament properties using a variety 
of statistical or topological models of large scale surveys. 
These procedures have estimated the width of filaments 
to be 1-2 Mpc (Ratcliffe et al 1996; Doroshkevich et al 
2004) and tens of Mpc in length (Bharadwaj, Bhavsar & 



Sheth 2004). 

With these limited direct constraints, TV-body simu- 
lations have played an instrumental role in our under- 
standing of the growth of the cosmic web. Dolag et al 
(2006) studied such a simulation and found that many of 
the filaments have density profiles that scale as r~ 2 -r~ 3 
with temperatures between ~10 5 -10 7 K. Similar results 
are found by Colberg et al (2005). Harford et al (2011) 
studied a similar filament at z ~ 5 and found slightly 
lower temperature 1 — 2 x 10 4 K, but with similar central 
densities of a few 10-1000 times the average background 
density. 

While star formation in a Milky Way size galaxy is rel- 
atively constant it is a different matter for dwarf galax- 
ies. The star formation per unit mass in these galaxies 
is just as intense as that seen in larger spiral galaxies 
without the expected dynamical causes, such as spiral 
arms and bars {e.g. Vanzi et al 2009). The star forma- 
tion is episodic and driven by mechanisms that are poorly 
understood, although mergers, tidal interactions, and ac- 
cretion from filaments are possible causes {e.g. Cumming 
et al 2008; James et al 2010). In addition, the metallicity 
of these galaxies is lower relative to larger star forming 
galaxies and they may provide insights to star formation 
during earlier formation epochs (Meier et al 2002). 

NGC 5253 represents a prime example of such a galaxy. 
It is a companion galaxy to M83 and is nearby at a 
distance of 3.8 Mpc (Gibson et al 2000). Recent ob- 
servations have uncovered multiple "super-star clusters" , 
which are very bright with Ljr ~ 10 9 L© {e.g. Gorjian, 
Turner, & Beck 2001) and M v ~ -10 to -11 (Meurer 
et al 1995; Gorjian 1996). In addition, radio observa- 
tions of CO(2-l) emission coincident with a dust lane 
and with properties consistent with inflow (Meier et al 
2002). However, the nature of the starburst and how 
large a part the infall plays remains an open question. 

It is with these observations in mind that we present 
the results from a series of 2D numerical simulations that 
explore the evolution of an idealized intergalactic fila- 
ment. We concentrate on the chemical and molecular 
makeup of the filament as it evolves from an initially 
atomic ionized state. At high temperatures, the cooling 
is dominated by atomic species. At lower temperatures, 
on the other hand, cooling from molecules and metals is 
more important. Here we track the formation and evo- 
lution of the key moleuclar species H2, CO, OH, and 
H2O, that contribute most to the cooling at low temper- 
atures. These results then allow us to asses the ability 
for filaments to cool to low temperatures and fragment 
into dense clouds. Finally, we extend our simulations 
to model a filament near a dwarf galaxy, comparable to 
NGC 5253. 

The structure of this paper is as follows. In §2 we in- 
troduce our chemistry and cooling routines and discuss 
the implementation and associated tests. In §3 we out- 
line our filament model and in §4 we present our results 
for isolated filaments. In §5 we present results for fila- 
ments forming near a dwarf galaxy, and a discussion and 
conclusions are given in §6. 

2. NUMERICAL METHOD 

All simulations were performed with FLASH version 
3.3, a publicly- available multidimensional adaptive mesh 
refinement hydrodynamic code (Fryxell et al 2000) that 



solves the Riemann problem on a Cartesian grid using 
a directionally-split piecewise parabolic method (PPM) 
(Colella & Woodward 1984; Colella & Glaz 1985; Fryx- 
ell 1990). Two new capabilities needed to be added 
to the public version of FLASH to properly simulate 
our filaments: a non-equilibrium chemistry package that 
tracks the formation and evolution of important molecu- 
lar coolants, and a cooling routine that takes into account 
the cooling from each of these molecules. In this section 
we describe the numerical implementation and tests for 
each of these new capabilities. 

2.1. Chemistry 

Even when only a few elements are considered, the gas 
phase chemistry of enriched gas can be quite complex, 
with hundreds of species and thousands of individual re- 
actions {e.g. Le Teuff et al 2010; Semenov et al 2010). 
While such large networks are important for many stud- 
ies of the interstellar medium (ISM), they are imprac- 
tical when coupled with large hydrodynamical simula- 
tions. We therefore used a substantially smaller network 
that tracks and evolves the atomic and molecular species 
made up of the elements that are most important in the 
thermodynamics of our system: hydrogen, helium, car- 
bon, and oxygen. 

2.1.1. Implementation 

In particular, we have implemented the chemical net- 
work presented in Glover et al (2010), which tracks the 
impact of 218 separate chemical reactions on 32 species: 
atomic hydrogen (H, H + , H _ ), atomic helium (He, He + ), 
atomic carbon (C, C + , C _ ), atomic oxygen (O, + , 
O - ), molecular hydrogen (H2, H^ , H^"), molecular car- 
bon (C2), molecular oxygen (O2, O^ ), key molecules con- 
taining combinations of these elements (OH, OH + , CO, 
CO+, CH, CH+, CH 2 , CH+,CH+, HCO+, HOC+, H 2 0, 
H20 + , H 3 + ), and free electrons (e~). As the binding 
and ionization energies of each species are important for 
the overall energy budget of the gas, we assigned each 
species with a energy, E v = B.E. — I.E, which is sim- 
ply the ionization potential subtracted from the binding 
energy, as summarized in Table [1] 

In addition to gas phase collision reactions, our net- 
work included the effect of an ambient ultraviolet (UV) 
radiation field. The photochemical rates assume a stan- 
dard interstellar radiation field from Draine (1978) which 
has a field strength of Go = 1.7 in Habing (1968) units. 
This corresponds to 6.5xl0~ 4 FUV photons per cm -3 . 
We introduced a coefficient, J21 = G/Go that allowed 
us to alter this background field, assuming that the field 
scales linearly with Go- For example, a J21 =0 removes 
the background field while a value of J21 = 1.0 gives the 
standard Habing field. 

Here we briefly describe our implementation of this 
network. First, we labele each species by an index i, 
such that species i had Z, protons, Ai nucleons, and a 
mass density of pi. From this we define a mass fraction 
as Xi = pi/ p, where p = J^i Pi i s the overall gas den- 
sity, and define the molar mass fraction as Yi = Xi/A; b . 
A continuity equation for each species is then given as 
Yi = Ri, where Ri is the total reaction rate due to all 
reactions. We then follow the evolution of every species 



TABLE 1 

Summary of species data. The first column gives the 

species name, the second column gives the ionization 

potential, third gives the dissociation energy, fourth 

gives e„, and the last column gives the ratio of specific 

heats. all energies are given in units of ev. those species 

denoted with a single asterisk do not have published 

dissociation energies. finally, hoc 4 ", which we denote 

with the double asterisk, does not have any published 

ionization or dissociation energies, and we assume they are 

the same as for hco+. 



Species 


LP (cV) 


D.E (cV) 


E„ (cV) 


7 


H 


0.00 


0.00 


0.00 


1.66 


H+ 


13.60 


0.00 


-13.60 


1.66 


H" 


0.00 


0.77 


0.77 


1.66 


He 


0.00 


0.00 


0.00 


1.66 


He+ 


24.60 


0.00 


-24.60 


1.66 


C 


0.00 


0.00 


0.00 


1.66 


C+ 


11.27 


0.00 


-11.27 


1.66 


C" 


0.00 


1.26 


1.26 


1.66 





0.00 


0.00 


0.00 


1.66 


o+ 


13.60 


0.00 


-13.60 


1.66 


o- 


0.00 


1.46 


1.46 


1.66 


H 2 


0.00 


4.48 


4.48 


1.40 


H+ 


15.43 


4.48 


-10.95 


1.40 


16.30 


6.30 


-10.00 


1.31 


c 2 


0.00 


6.21 


6.21 


1.40 


o 2 


0.00 


5.12 


5.12 


1.40 


0+ 


12.06 


6.66 


-5.40 


1.40 


OH 


0.00 


4.39 


4.39 


1.40 


OH+ 


13.00 


5.10 


-7.90 


1.40 


CO 


0.00 


11.90 


11.90 


1.40 


CO+ 


14.01 


8.34 


-5.67 


1.40 


CH 


0.00 


3.47 


3.47 


1.40 


CH+ 


10.64 


0.00 


-10.64 


1.40 * 


CH 2 


0.00 


4.00 


4.00 


1.31 


CH+ 
CHJ- 


10.40 


4.00 


-6.40 


1.31 


9.83 


0.00 


-9.83 


1.31 * 


HCO+ 


9.88 


0.00 


-9.88 


1.31 * 


HOC+ 


9.88 


0.00 


-9.88 


1 Q1 ** 


H 2 


0.00 


0.00 


0.00 


1.31 


H 2 0+ 


12.61 


5.11 


-7.50 


1.31 


H s O+ 


0.00 


0.00 


0.00 


1.31 * 


c~ 


0.00 


0.00 


0.00 


1.66 



in our simulation except for free electrons which we di- 
rectly calculate using conservation of the total electron 
density. 

Due to the complex ways that chemical reaction rates 
depend on temperature and the wide range of abun- 
dances possible, the resulting rate equations arc 'stiff'. 
This requires implicit or semi-implicit methods to prop- 
erly solve these equations. To that end, we have imple- 
mented a fourth-order accurate Kaps-Rentrop, or Rosen- 
brock method (Kaps & Rentrop 1979), as described in 
Gray & Scannapieco (2010; hereafter GS10). As the 
species evolve, the temperature will change as the inter- 
nal energy is altered from recombinations, ionizations, 
and dissociations. Since the reaction rates are strong 
functions of temperature, the network can become un- 
stable if too large of a time step is used. To ensure the 
stability of the network while running the simulation at 
the hydrodynamic time step, we use the same subcycling 
method from GS10. A chemical time step was calculated 
for each species as 



7"chcm,'. 



^chem" 



Yi + 0.1Y H+ 
Yi 



(1) 



where a c hcm is a constant determined at runtime and 
controls the maximum abundance change allowed, which 
we default to 0.1. We used the current value for the 
species abundances and the change in abundances (Yi) 
were calculated from the ordinary differential equations 
at the current temperature. Finally, we included a small 
fraction of Y H + as a buffer to prevent species with very 
small abundances, but changing rapidly, from causing 
prohibitively small time scales and therefore excessive 
subcycling. 

At the beginning of the chemistry cycle, the smallest 
species time scale was compared to the hydrodyanmic 
time step. If the hydro time is smaller than the chem- 
istry time step, no additional subcycling is done and the 
network was advanced as normal. If the chemistry time 
step is smaller than the hydro time step, then the net- 
work is advanced using the chemistry time step. The 
newly updated species and temperature are used to re- 
calculate the chemistry time step, which was then com- 
pared to the remaining hydro time. This continues until 
the species have advanced a full hydro time step. 

Finally, there are two special cases where the chemistry 
can be simplified. At temperatures greater than 1.7x 10 7 
K, we assumed that all species are atomic and ionized 
and no additional chemistry is done. At temperatures 
greater than l.OxlO 4 K but less than 1.7xl0 7 K, we 
assume that all molecules are dissociated and added to 
the atomic neutral species (e.g. Xh = Xh + Xh 2 + ...). 
A chemistry step is then performed. 

2.1.2. Chemistry Test 

As a test of our implementation of these coupled equa- 
tions, we used a procedure similar to that presented in 
GS10 and Glover et al (2010). In order to have a com- 
mensurate dataset to compare against, we use reaction 
rates given in the UMIST RATE06 database (Woodall et 
al 2007), and restricted our attention to two-body inter- 
actions in the absence of a background UV dissociating 
field. Out of the 4604 reactions in RATE06, 134 were 
used that matched the reactions given in Glover et al 
(2010). 

We completed this complementary network by adding 
rates that appear in Glover et al (2010) and not in 
RATE06. In total, this gave a total of 144 reactions. In 
particular we added reactions for the positive-negative 
ion recombination of O - and H + (Rxn 108 in Glover 
et al 2010), the formation of H 3 0+ from H 2 and H+ 
(Rxn 141), the formation of Hj from H~ and H + (Rxn 
15), the formation of H^~ from H and H + (Rxn 4), the 
charge exchange reaction between H2 and H + (Rxn 7), 
ion molecule reactions between H30 + and H (Rxn 55), 
between + and H2 (Rxn 69), and between CO and Hc + 
(Rxn 105), and the electron attachment of HCO + (Rxn 
130). 

Next, we compare the reaction rates presented in 
Glover et al (2010) to those in RATE06. When these 
rates are appreciably different, we alter those rates to 
match those found in RATE06. Specifically, we changed 
the dissociation of H2 by H (Rxn 9) and of H2 by H2 (Rxn 
10), the radiative recombination of H + (Rxn 12) and He + 
(Rxn 17), the formation of CO by CH and O (Rnx 38), 
the formation of O2 by OH and O (Rxn 47), the forma- 
tion of H^ by H^ and H2 (Rxn 54), the dissociation of 



CO by He + (Rxn 104), the dissociative recombination 
of Hj (Rxn 110-112), the dissociative recombination of 
H^ (Rxn 120-122), and the dissociative recombination of 
H 3 0+ (Rxn 123-126). 

Having matched our FLASH network with a restricted 
set of RATE06 reactions, we then ran a series of simple 
tests. We use a series of single zone models to com- 
pare the chemical evolution over a wide range of intital 
temperatures and densities. In each model, the ther- 
modynamic variables are kept constant. To advance the 
RATE06 reaction network, we use the open source chem- 
istry code AstrocheniJ. Models were run for a total of 
10 8 years with hydrogen number densities ranging be- 
tween 10" 2 to 10 6 cm -3 and temperatures ranging be- 
tween 10 2 to 10 4 K. In all cases the initial abundances of 
each species relative to the total hydrogen density were: 
n Hc = 0.08, n c += 8.2x 10" 5 , n + = 1.5x 10~ 5 , n m = 
0.99, and n-n = 0.01. We then ran similar models with 
our modified FLASH code, starting each model with a 
small initial time step (~ 10 5 s) which was allowed to 
increase to the hydrodynamic time step. 

Figure [1] shows a comparison between the derived 
chemistry network evolved with Astrochem and the 
FLASH network where the temperature is constant at 
100 K and a hydrogen number density of ijh = 10~ 2 
cm~ 3 . While the abundances change over many orders 
of magnitude through the evolution time, our FLASH 
always closely match the Astrochem network. Similar 
results are found for the range of conditions considered. 

Another simple test was to ensure that the chemi- 
cal abundances scaled appropriately with density. The 
timescale for formation of a given chemical species can be 
estimated as T c h C m ~ n/h, where n and n are the instan- 
taneous species number density and its rate of change, 
respectively. At low densities where three body reac- 
tions are unimportant, fi is given by kn 2 , where k is the 
formation or destruction rate, a function of only temper- 
ature. This gives a chemistry timescale that is inversely 
proportional to the overall density. We ran several tests 
in which where we compared the time evolution of our 
chemical species at a constant temperature and changed 
the initial mass density. Although not shown here, the 
species always evolved as expected. 

2.2. Cooling 

As the species evolve, the temperature of the gas 
changes due to recombinations, ionizations, and disso- 
ciations, as described above. In addition, at low tem- 
peratures (ilO 4 K), the cooling from molecular species 
becomes important. Since both the reaction and cooling 
rates are temperature dependent, it is important for both 
of these processes to remain coupled. Therefore, once a 
chemical subcycle was performed, the cooling routines 
were immediately called. This ensured the stability of 
both the chemistry and cooling routines. 

At temperatures >10 4 K, the majority of radiative 
cooling comes from atomic species, with bremmstrahlung 
radiation becoming important at temperatures >10 7 K. 
Below 10 4 K, many cooling channels become important, 
primarily from atomic metal lines from carbon and oxy- 
gen, as well as molecular line cooling from H2, CO, H2O, 

3 http://smaret.github.com/astrochem/ 
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Fig. 1. — Comparison of chemical abundances between As- 
trochem and FLASH. In each panel the lines show the abundances 
from FLASH while the circles are from Astrochem. The y-axis 
of each panel is the logarithm of time in units of years and the 
number density of the given species divided by the total number 
density of hydrogen. Top Right: The dashed (red) lines show H + , 
the solid (blue) lines show H, the fine-dashed (magenta) lines show 
He+, and the dotted (green) lines show He. Top Left: The dashed 
(red) lines show O, the solid (blue) lines show C, the fine-dashed 
(magenta) lines show + , and the dotted (green) lines show C + . 
Bottom Right: The dashed (red) lines show H2, the solid (blue) 
lines show CO, and the dotted (green) lines show OH. Bottom Left: 
The dashed (red) lines show CH, the solid (blue) lines show CH2, 
the dotted (green) lines show O2, and the fine-dashed (magenta) 
lines show C2- In general, the solid lines and circles always overlap 
each other and it is very hard to differentiate between them. 

and OH. In particular, carbon monoxide (CO) is ex- 
tremely important at temperatures below 200 K where 
molecular hydrogen (H2) becomes inefficient (Smith et 
al 2008). Finally as in Glover et al (2010), we included 
heating from the photoelectric effect, H2 UV pumping. 
In this case, the total cooling function becomes, 



A 



Total 



A 



H,Ho 



A 



Molecules 



A 



Metals 



A 



Heating 



(2) 



such that the total cooling function in units of ergs s _1 
cm~ 3 , is simply the summation of the cooling rate from 
each cooling channel minus sources of gas heating. The 
introduction of carbon and oxygen into the chemistry 
network allows for the creation of many molecules that 
are important for the thermal evolution of the gas. As we 
have expanded the chemistry network, we also expanded 
the cooling routine presented in GS10 to include these 
important coolants. Much of this cooling follows from 
Glover et al (2010), which we summarize below. 

2.2.1. Atomic Cooling 

The treatment of the high temperature hydrogen and 
helium cooling remained unchanged from GS10. The 
cooling rate is calculated using CLOUDY (Ferland et al 
1998) assuming only collisional ionization in the Case B 
limit. 



To account for metal-line cooling, we used the tabu- 
lated results from Weirsma et al (2008), assuming stan- 
dard solar rations scaled by the local metallicity. These 
rates assume local thermodynamic equilibrium and are 
defined from 10 2 to 10 9 K. Note that while cooling by 
atomic hydrogen and helium are important above 10 4 K, 
below this temperature atomic cooling comes primarily 
from neutral oxygen and neutral and singly ionized car- 
bon. 

2.2.2. Molecular Cooling 

Hi cooling — Much of the H2 cooling rates are taken 
from Glover et al (2008) and again remain unchanged 
from GS10. Each cooling rate has the form: 

A i:j = runjX^j, (3) 

where A^- is the energy loss per volume due to species i 
and j, rii and rij are the number densities of each species, 
and Xij is the cooling rate in ergs cm -3 s _1 . Cooling 
rates for the collisional excitation between H2 and H, 
H2, H + , and e~ and between H^ and H or e~, are all 
included. 



CO and H2O cooling — To implement CO and H2O cool- 
ing, we followed Glover et al (2010). Here the au- 
thors used the tabulated rotational cooling functions 
from Neufeld & Kaufman (1993) and Neufeld, Lepp, & 
Melnick (1995). As above, the cooling rate is defined 
as a rate coefficient L, with units of erg cm 3 s _1 such 
that the cooling rate per unit volume per unit time is 
A=Ln(H2)n(M), where n(i? 2 ) is the number density of 
H2 and n(M) is the number density of the other molec- 
ular species. 

The cooling rate coefficient is given as a function of 
the H2 number density, the kinetic temperature, and a 
parameter that relates the species number density to the 
local velocity field N(m), defined as: 



N(m) = 



n(m) 



(4) 



The rate coefficients for CO and H2O cooling are com- 
puted as a four-parameter analytic fit of the form: 



1 
I 



1 

T 



n(H a ) 



^LTE 



1 

T 



n(H 2 ) 

%/2 



1 



"1/2-^0 

£lte 



(5) 



Here Lq represents the cooling in the low density limit, 
Lute is the cooling rate when the rotational level popula- 
tions are in local thermodynamic equilibrium, and TI1/2 is 
the H2 number density when L — O.BLq. Lq is a function 
of temperature only while the other three parameters 
(Llte, ^i/2) an d a) are functions of temperature and 
the effective column density per unit velocity as defined 
in eq. (g}. 

Rotational cooling for H2O the cooling rate coefficient 
was tabulated between 10 K < T <, 4000 K and for col- 
umn densities between 10 < log 10 (A r (H 2 O)) < 19, where 

iV(H 2 0) in units of cm -2 per km s -1 . We assumed a 
constant ortho-para ratio of 3:1 as Neufeld et al (1995) 
who define cooling parameters for both states of H 2 0. 
Similarly, the rotation cooling of CO is defined for tem- 
peratures between 10 K < T < 2000 K and for column 
densities between 14.5 < log 10 (iV(CO)) < 19. 



We also introduced an artificial temperature floor at 
10 K and turn off any radiative cooling at these tem- 
peratures. To properly model gas to temperatures below 
this value, one would be required to not only extend the 
results from Neufeld et al (1993; 1995) but also include 
other physical processes, such as the CO freeze-out (e.g. 
Lee et al 2004) or, if included in the chemical network, 
the decrease in dust temperature as the extinction in- 
creases (e.g. Goldsmith 2001). 

For gas with temperatures above the tabulated values, 
we simply adopted cooling rates that correspond to the 
highest tabulated temperature. If the column density 
was within the tabulated range, then we linearly inter- 
polated along column density to determine our cooling 
rate. Although at such high temperatures, it is expected 
that these coolants are quickly dissociated. 

A similar process was used when exceeding the tabu- 
lated column density. If a zone was denser than tabu- 
lated, we used the highest values given for the cooling- 
rate and, if possible, interpolated along temperature. If 
a zone was more rarified than tabulated, then we used 
the lowest tabulated values. At very high densities, this 
procedure may overestimate the expected cooling some- 
what. Although, in our simulations the abundance of 
these coolants was negligible at such low densities, and 
they did not contribute much to the overall cooling. At 
lower densities however, the cooling rates naturally ap- 
proach the optically thin limit and it is unlikely to intro- 
duce much error in our simulations. 

Finally, in their calculation of the cooling functions, 
Neufeld et al (1993; 1995) assumed that only collisions 
between the target coolant and H2 are important. How- 
ever, as pointed out by Glover et al (2010), if the gas 
is not completely molecular then the collisions between 
atomic hydrogen and electrons also contribute. To ac- 
count for this, we replaced «h 2 hi eq. [5]and in the total 
cooling rate, with an effective number density, n e g. For 
CO rotational cooling the effective number density is 



™eff,CO,rot ="H 



M^ 



n H 

1.3 x 10" 8 cm 3 s- 1 
an 2 v e 



(6) 



where hh 2 , uh, and n c are the number densities of 
H 2 , hydrogen, and electrons respectively, cth = 2.3 x 



10 



-15 



cm 2 , ct H2 = S.SxlO-^T/lOOOK)- 1 ^ 



and 



w e =1.03xl0 4 y/T(K) cm s" 1 . For H 2 rotational cool- 
ing the effective number density is 



™cff,H 2 0,rot — "-H 2 



KH 2 



(7) 



where k e = dex[-8.020 + 15.740/T 1 / 6 - 47.137/T 1 / 3 + 
76.648/T 1 / 2 - 60.191/T 2 / 3 ] and k H2 = 7.4xl0~ 12 T 1 ' 2 
cm 3 s _1 . These equations are taken from Meijerink & 
Spaans (2005) while the formula for k e is taken from 
Faure, Gorfinkiel & Tennyson (2004). 

Vibrational cooling from CO and H2O , which is impor- 
tant at high temperatures and densities, is also presented 
in Neufeld et al (1993). These authors provide a simpler 



two parameter fit of the form 



1 
L 



1 



UH a 



'LTE 



(8) 



Analytical functions are provided for Lq for CO and H2O 
as well as tabulated values for Llte- Vibrational cooling 
rates for both coolants are tabulated for temperatures 
100 K <T < 4,000 K and for column densities of 13 

< log 10 (JV) < 20. As above, «h 2 is replaced with up- 
dated values from Meijerink & Spaans (2005). For CO 
vibrational cooling, the effective number density is 



"■cff.co.vib = nn 2 



50«h 



L 



CO,c 
^CO,0 



(9) 



where L co , c = 1-03 x lO" 10 (^^'^exp ( 



-3080 > 



L 



CO,0 



1.14 x 10" 14 exp (^y#) exp 1 



- — ), and 
3S§2). Similarly 



for H2O vibrational cooling, the effective number density 
is 



«eff,H 2 o,vib = "h 2 + 10n H + I , 2 ,e I II, 



L 



H 2 O,0 



where Lh 2 o,c = 2.6 x 
L H ,on = 0.64 x 10~ 14 



^H 2 O,0 



exp 



10 -6 T -l/2 

- 47 - 5 ^exp 



-2325 i 



(10) 



and 



exp ^ — jr- 

2325 ^ For tern- 



jTJ3~ J i^}/ V T 

peratures and effective column densities outside of the 
tabulated range, we calculated cooling rates as we did 
for the rotational case. 

OH cooling — To model the cooling from OH we used 
the tabulated results from Omukai et al (2010), who 
computed the OH cooling rate in the same manner as 
Neufeld et al (1993; 1995) and presents values to fit eq. 
([5]). The rate is valid in a temperature range between 30 
K < T <J300 K and for effective column densities of 10 
< logio(N (OHj) < 18. For temperatures and densities 
outside this range, we calculated the appropriate cool- 
ing as above. Here we assume that only direct collisions 
between OH and H2 determine the cooling rate. 

Heating — To complete our treatment of the thermal evo- 
lution of the gas, we included heating by the photoelec- 
tric effect and H2 UV pumping as modeled by Glover et 
al (2010). We also included a rate that takes into ac- 
count the photoionization of hydrogen. When hydrogen 
is collisionally ionized, energy is taken from the gas to 
overcome the binding energy of the atom and lowering 
the internal energy of the gas. However, if the atom is 
photoionized, it does not alter the internal energy. We 
model this as a heating term of the form: 



r H + =Ar H +£H 



t H+7->H+ 



£"fcj + fcH+7^H+ 



(11) 



where AYfj+ is the difference of H + is a time step, £h is 
the binding energy of hydrogen, fcn+ 7 ^H+ is the rate at 
which hydrogen is photoionized, and S"fej is the sum- 
mation of all the collisional ionization rates that form 
H + . Similar equations can be written for other atomic 
species, such as helium and carbon, however, hydrogen 
always provides the dominant heating term. 



2.2.3. Cooling Tests 

Having tested the chemistry package, we turned our 
attention to the cooling routines. Several of the cooling 
channels have been tested and used in previous appli- 
cations, such as the high temperature atomic hydrogen 
and helium cooling and molecular hydrogen cooling (e.g. 
GS10). For cooling from CO, OH, and H2O, we simply 
ensured that the interpolation between points in the ta- 
ble was smooth and returned the table values under the 
correct conditions. 

3. MODEL FRAMEWORK 

Once we had implemented the required chemistry and 
cooling routines into FLASH, we were then able to turn 
our attention to the structure and composition of the 
filament. We assume a Cold Dark Matter (CDM) cos- 
mology with parameters, fl\ = 0.7, f^o = 0.3, fib = 0.045, 
and h = 0.70, (e.g. Spergel et al 2007), where h is the 
Hubble constant with units of 100 km s _1 Mpc" 1 , ft\, 
J7o, and fib are the vacuum, total matter, and baryonic 
matter densities in units of the critical density. For our 
value of h, the critical density p cr u is 9.30xl0~ 30 gm 
cm~ 3 , which corresponds to a mean baryonic density of 
p b =4.18xl0~ 31 gm cm~ 3 . 

3.1. Filament Model 

Building on previous analytic and numerical results we 
model the filament as a combination of dark and bary- 
onic matter with a density profile (Stodolkiewicz 1963; 
Ostriker 1964; Harford et al 2011) of 



P(i = r/r s ) 



(1 + tO 2 SmCm 



(12) 



where p c is the central gas density and r s is the scale 
radius. If the filament were made of pure gas r s it would 
be in hydrostatic equilibrium if 

1/2 



( k b T 



(AttGpc) 1 / 2 \4irGp c pm p 



(13) 



where c s is the sound speed, kb is Boltzmann's constant, 
p is the average molecular weight, m p is the mass of a 
proton, and G is the gravitational constant. However, 
to account for the presence of dark matter we need to 
redefine this quantity, making use of the equation of hy- 
drostatic equilibrium 



V$ 



1 



Ptotal 



-VP 



(1 + 0)P, 



-Vp«, 



(14) 



where P is the thermal gas pressure, p g is the baryonic 
gas density, and j3 = (fig — Oj,)/f2(,, where the constants 
are defined above. Using this equation we can also define 
an effective sound speed as 



c 2 je// ^c 2 /(l + /3). 
Combining this with Poisson's equation 

V 2 $ = -47rGptatal = -4ttG(1 + 0)p g 
alters the definition of the scale radius to 



(15) 
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kpc, (18) 



where we re-express the central density (/9 S , C ) as Apb 
where A is a constant defined at runtime that represents 
the overdensity of the filament. 

The composition of the filament is assumed to be ini- 
tially ionized and atomic with abundances of X^ = 
0.76, A Hc + = 0.24, X c+ = 1.693xlO- 3 Z frac , and X Q+ 
= 5.056 xl0~ 3 Zf rac , where Zf lac is the metallicity of the 
gas in solar units and X[ is the mass fraction of the in- 
dividual species. 

The long standing discrepancy between models of 
chemical evolution in galaxies and observations has given 
rise to the well known G Dwarf problem, first noticed by 
van den Bergh (1962) and Schmidt (1963). Simple chem- 
ical evolution models produce a population of metal-poor 
stars that are not seen in observations. One solution is to 
assume that the gas has some initial metallicity, perhaps 
provided by outflows from dwarf galaxies at high redshift 
(e.g. Scannapicco 2005). To that end, we assumed that 
in these low redshift cases that the filament gas has been 
enriched to Z = 0.1 Z Q . 

3.2. Gravity and the Jeans Limit 

To account for the gravitational potential due to the 
surrounding dark matter, we calculate the gravitational 
acceleration of the gas as, atotai = asc + a DM, where asa 
is the self gravity of the gas as computed by the Poisson 
solver within FLASH (Rickcr 2008), and the acceleration 
from the dark matter is given as 



odm = -V$ = -2c z s f3 Vln I 1 + ^ 



(19) 



At the beginning of each simulation, the baryonic mat- 
ter was assumed to be in hydrostatic equilibrium with 
the dark matter. A test of this gravity scheme showed 
that the filament remains in hydrostatic balance for many 
free-fall times in the absence of any chemistry or cooling. 
A common problem that arises in hydrodynamical sim- 
ulations is artificial fragmentation of perturbations that 
comes from purely numerical sources. Truelove et al 
(1997) studied this problem and found that to avoid this 
numerical artifact, the spatial resolution of the simula- 
tion must be smaller than the Jeans length, defined as 
(Jeans 1902), 

Aj = — TT7*. (20) 



(G P ) 



1/2' 



where c s is the local sound speed. Truelove suggests that 
the ratio between the Jeans length and the spatial resolu- 
tion, known as the Jeans number (Nj = Xj/Ax), should 
be greater than 4 for hydrodynamical simulations, al- 
though, cosmological simulations have found that a value 
of 7 is required (Ceverino et al 2010). 

With the advent of adaptive mesh refinement codes, 
this requirement becomes easier to manage since one can 
refine only those regions that require it. However, this 
is feasible only to a point. When the Jeans length gets 
below a few times the resolution of even the highest AMR 



level, we introduce artificial pressure support of the form, 
P art = C 2 Gp 2 (Ai min ) 2 //«, (21) 

where C = 10, is a dimensionless constant proportional 
to the number of cells being smoothed over, p is the den- 
sity, and Ax m i n is the minimum cell size (Machacek et 
al 2001; Agcrtz et al 2009; Ceverino et al 2010; Klar 
& Miicket 2012) which is used in place of the thermal 
pressure in further hydrodynamic evolution. 

Naturally, once this condition is met, we can no longer 
follow the small scale density evolution of that region 
since the artificial pressure will act to smooth out these 
perturbations. Therefore, in the analysis that follows, we 
will describe these regions in terms of average quantities 
instead of specific values. 

Finally, the use of a pressure floor would cause the 
temperature to rise to unphysical values if not accounted 
for in the equation of state. To take this into account, we 
tracked the fraction of the pressure that comes from the 
artificial source and corrected the temperature for use in 
the chemistry and cooling routines. 

It should also be noted that we begin each simulation in 
hydrostatic equilibrium, but not in thermal equilibrium. 
The reason for this is two-fold. First, the cooling expe- 
rienced from hydrogen-helium recombination and metal- 
line cooling is very strong, and second, the possible heat- 
ing terms are very small. 

The number of possible heating mechanisms is very 
limited and largely depend on the UV or cosmic ray back- 
ground. In the absence of these conditions, the heating 
is negligible. Even with the presence of a nominal back- 
ground, the photoelectric heating is only marginal with a 
rate of w l.Ox 10 _26 n ergs s _1 cm~ 3 , for a 10 6 K gas, an 
efficiency of w 1, and G — 0.01, where n is the number 
density. Similarly for cosmic ray heating, the nominal 
value for the cosmic ray ionization rate is w l.OxlO -17 
s _1 (Dalgarno 2006), which leads to a heating rate of 
s=s 10~ 28 n ergs s _1 cm~ 3 . Both of these rates are much 
smaller than the cooling mentioned above. 

4. RESULTS 

Each of our simulations is carried out in two dimen- 
sions with cylindrical coordinates with the r and z axis 
spanning [0,5] r s . If the ratio of the first and second 
derivatives of density, temperature, or pressure exceeded 
0.6 divided by the cell size, then the block was marked 
for refinement. Conversely, if the value of this ratio for 
all three variables was less than 0.20 over the cell size, 
it was marked for dcrcfincment. If a filament were to 
lose all thermal pressure support, it would collapse on a 
free-fall timescale given by 



Tff 



3vr 



32G/3A Pc 



0.9 Gyrs 



(22) 



and, in addition to the criteria above, we forced derefinc- 
ment in regions that were greater than 2.0 r s from the 
center of the cylinder at times greater than 0.5 Tff. Each 
simulation was run with a maximum of 11 refinement lev- 
els, which gave a large dynamic range of 8192 and a min- 
imum cell size of 13.2 pc in our fiducial case. We imposed 
a reflecting boundary condition along the central axis and 
a diode boundary condition on the right, which allowed 
gas to leave the simulation area and prevents inflow from 



TABLE 2 
Summary of completed simulations. 



Name 



Ovcrdcnsity Temp. (K) R 3 (kpc) UV Background Metallicity (Zq) Tg (Gyrs) 



Fiducial 


(FID) 


2000 


1.0x10 s 


21.6 


Ldcns 




200 


1.0x10 s 


68.3 


HightT 




2000 


l.OxlO 7 


68.3 


UV1 




2000 


1.0x10 s 


21.6 


UV2 




2000 


1.0x10 s 


21.6 


Highz 




500 


1.2xl0 4 


0.83 



0.0 
0.0 
0.0 
1.0 
0.1 
0.0 



l.OxlO- 1 
LOxlO" 1 
LOxlO" 1 
l.OxlO" 1 
LOxlO" 1 
LOxlO" 3 



0.895 
2.830 
0.895 
0.895 
0.895 
0.385 



these boundaries. We also impose diode boundary con- 
ditions at the top and bottom of the simulation domain. 
A summary of the simulations performed is given in 
Table [2] For each run we give the initial temperature, 
overdensity of the cylinder, corresponding scale radius, 
strength of the UV background, initial metallicity of the 
gas, and associated free-fall time. 

4.1. Hydrodynamic Evolution 

The initial density, pressure, and temperature profiles 
for our fiducial case are shown in Figure [3J The cool- 
ing experienced by this gas largely determines its hy- 
drodynamic evolution. Initially, cooling is dominated by 
atomic cooling from metals and hydrogen and helium. 
The atomic cooling time scale, i.e. the time required to 
completely cool the gas from an initial temperature, T, 
is given by 



l.bnkhT 



7~cool 



n H n e (A H + Z iTac A M ) 



00, 



(23) 



where n, n#, and n e are the total number density, hy- 
drogen number density, and electron number density re- 
spectively, and Ah and Am are the atomic and metal- 
line cooling rate respectively. For our initial tempera- 
ture of 10 6 K, this time scale is ps 1.8 Gyrs. However, 
cooling becomes more efficient at slightly lower temper- 
atures ( ~ 10 5 K, e.g. Sutherland & Dopita, 1993) until 
the hydrogen-helium cooling floor is reached at 10 4 K, at 
which point the timescale is a fraction of a Gyr. 

In Figure. [3l we show ID profiles of the evolution of our 
fiducial filament. The top left panel shows the density 
profile, the top right shows the temperature profile, the 
bottom left shows the H2 mass fraction, and the bottom 
right panel shows the CO mass fraction. As the filament 
gas begins to cool, the initially ionized gas rapidly re- 
combines and cools to near the hydrogen-helium cooling 
floor. At this point, the gas begins to collapse toward 
the center of the filament which increases the density, 
and therefore, the cooling experienced at the center. The 
relative importance of each cooling term is shown in the 
top left panel of Figure HI Early on, cooling is dominated 
by atomic hydrogen-helium and metal-line cooling. How- 
ever, once the cloud becomes dense enough for molecules 
to form, their cooling channels become as important as 
the atomic cooling. 

As the gas cools, shocks are formed as the infalling gas 
interacts with this cooler interior gas. These create peaks 
in the density profiles and allow for an increase in the 
abundances in molecular species (cyan and magenta lines 
in Fig 13]). The shocks are then reflected off the central 
axis. However, these density peaks are quickly erased 
and the gas returns to a much more uniform distribution 



as shown by the dash-dotted and long-dashed lines in 
Figure [3 

The final state of the filament is shown by the long- 
dashed lines in Figure. [3] The gas has cooled and col- 
lapsed into a compact structure with a high central den- 
sity (p ps 1.0 x 10 -21 gm cm -3 ) and cold central temper- 
ature (T sa 10 K). In addition, the gas is enriched with 
molecular coolants with interior abundances of 10~ 3 for 
H2 and 10 -6 for CO respectively. These values should 
be taken as averages since at these temperatures we are 
forced to use our alternative pressure scheme to prevent 
unphysical collapse. 

4.2. Effect of Density 

The variable A controls the many of the properties 
of the filament, including its central density, size, and 
evolution timescale. More importantly, both chemistry 
and cooling rates are dependent on this density but with 
different scalings. As can be seen from eq. (|22l) . the free- 
fall timescale scales as p -1 / 2 while the cooling timescale 
(eq. |2"3")) goes as p _1 . Thus, in addition to our fiducial 
run, we ran a simulation where A was reduced by a factor 
of 10. 

Figure [5] compares the evolution of the density, temper- 
ature, and H2 mass fraction in this simulation with the 
fiducial run at similar evolutionary stages. In general, 
the low density and fiducial cases evolve in similar ways. 
Initially, the gas recombines and the filament begins to 
collapse, but the cooling in the lower density gas is not as 
efficient. This prevents shocks from forming and allows 
for a much more uniform collapse. Ultimately the out- 
come is the same, however, a dense core is formed at the 
center of the filament, with a large molecular fraction. 

The top right panel of Figure 3] shows the relative im- 
portance of each cooling term. Again, there is a strong 
similarity with the fiducial case. In each case, atomic 
hydrogen-helium and metal-line cooling dominates early 
on and, after an abundance of molecular coolants are 
formed, they become an important cooling channel. An 
important note is that while in the fiducial case the whole 
filament collapses within a single free-fall time (0.905 rg 
in the fiducial case), the low density case takes much 
longer (2.63 rjf), due to the increase in the cooling time, 
which scales as one over the number density. This can 
been seen when comparing the results in the top row of 
Fig|U In the fiducial case, for example, molecular cooling 
becomes important at w 0.7 Tg, while in the low density 
case, it does not become important until sa 2.3 rg . 

4.3. Effect of UV Background 

Next we explored the impact of a dissociating UV back- 
ground on the chemical and hydrodynamical evolution of 
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Initial simulation setup for our fiducial parameters. Top Left: Logarithmic density contours between 10 
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< p < 10" 



cm . Top Right: Initial logarithmic pressure contours between 10 < P < 10 ergs/cm . Bottom: ID profiles of the initial density 
(solid blue line) and pressure (dashed red line). 



the filament. We ran two simulations with different UV 
backgrounds, one with a full strength Habig field (UV1), 
denoted by a J21 value of 1.0 and a field with one tenth 
that strength (UV2). This field has the primary effect of 
dissociating molecular coolants after they have formed, in 
particular, molecular hydrogen (H2) and carbon monox- 
ide (CO). 

Figure [5] shows the ID profile comparisons between 
these runs. The first column shows our fiducial case, the 
second column shows UV1, and the third column shows 
UV2. Since the free fall time for these simulations is the 
same, each of the lines corresponds to the same evolu- 
tionary time. Comparing the density profiles between 
these run (first row of Figured]) shows that there is little 
difference. All three evolve with the same shock struc- 
ture developing and with similar final central densities. 



Slight differences can be seen in the temperature profiles 
where the fiducial case cools slightly faster than either 
UV1 or UV2 at late times. As expected, the greatest dif- 
ferences can be seen in the molecular hydrogen profiles. 
The stronger the UV background the lower the abun- 
dance of molecules. However, this has little impact of 
the final state of the filament, the cooling experienced, 
or even the final abundances found at the center of the 
filament. 

The middle row of Figure |4] shows the relative impor- 
tance of the various coolants for these runs. The primary 
difference is that while cooling from CO is important for 
FID, it is much less important for the UV runs. Cooling 
from OH, while marginal for FID is non-existent in UV1 
and UV2. The UV background also effects the time at 
which molecular cooling becomes important. The lower 
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Fig. 3. — ID profiles from the fiducial run. The x-axis in each panel is the logarithm of the radial position along the filament. Each 
y-axis is the logarithm of each variable. The top left panel shows the density profiles, the top right the temperature profiles, bottom left 
shows the H2 mass fraction, and the bottom right shows the CO mass fraction. Times are given in the legend in the density plot in units 
of the initial free fall time. 



the incident background, the earlier the cooling is impor- 
tant. 

4.4. Effect of the Initial Temperature 

The initial temperature of the filament not only de- 
termines the physical scale of the filament (eq. [l"3]) . but 
also the cooling rate experienced by the gas. To study 
the impact that the initial temperature has on the evolu- 
tion of the filament, we run a simulation with an initial 
temperature of 10 7 K (HighT). 

In fact, it is the cooling rate that completely deter- 
mines the fate of this filament. With the initial tem- 
perature so high, the cooling timescale becomes incredi- 
bly long. In fact, with an initial central number density 
of ^10~ 4 cm~ 3 and an initial temperature of 10 7 K, the 
cooling timescale is greater than 10 Gyrs, whereas the 
free-fall timescale is only 0.9 Gyrs. This allows the fila- 
ment to remain close to hydrostatic balance over several 
filament free fall times. 



4.5. High Redshift Filament Evolution 

Finally, we explored the evolution of a high-redshift 
filament. Recent simulations by Harford et al (2011) 
studied cosmological simulations at a z — 5 and found 
baryon-rich filaments that stretched between galaxies. 
They found that these filaments were well described by 
an isothermal cylinder with a typical overdensity above 
the cosmic mean density of 500, temperatures of ~ 1- 
2xl0 4 K, and nearly fully ionized gas. We use these 
values as initial conditions, also assuming a lower metal- 
licity of Z/Zq = 10~ 3 , in this case, since by this age in 
cosmic time, fewer stars have formed and have polluted 
the surrounding gas. 

Figure [7] shows the comparison between Highz and 
FID. As in the high temperature case, the filament does 
not collapse and instead relaxes into a hydrostatically 
balanced state. That initially ionized gas, again, quickly 
recombines which allows for the formation of a small 
abundance of molecular coolants. However, since the ini- 
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Fig. 4. — Summary of the density averaged cooling for each simulation. The top left panel shows the cooling terms for our fiducial run, 
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cooler, molecular cooling is more important much earlier. 
Finally, over much of the evolution of the filament, cool- 
ing from H2 is more important than metal-line cooling. 
However, this cooling is still not enough to allow the fila- 
ment to continue to continue its collapse after the initial 
atomic cooling stage. 



tial temperature is close to the atomic hydrogen-helium 
cooling limit, the increase in density of the filament dur- 
ing the atomic cooling stage is minimal. This means 
that densities are too low for molecular hydrogen, which 
is a much less efficient coolant, to take over at this stage 
and continue the collapse of the filament. Instead, the 
filament reaches a new quasi-steady state, which it will 
maintain for many dynamical times. 

The bottom right panel of Figure 2] the relative cooling 
experienced. This differs from FID in two aspects. First, 
since the metallicity is much lower in Highz than in FID, 
the relative contribution of metals is much lower. Sec- 
ond, since the Highz filament is initially much denser and 



5. A FILAMENT NEAR A DWARF GALAXY 

Dwarf galaxies are unique places in the study of star 
formation. Without the aid of normal dynamic causes, 
such as galactic bars or spiral arms (e.g. Vanzi et al 
2009), dwarf galaxies experience the same level of star 
formation per unit mass as their larger cousins. This 
intense, episodic star formation might be caused by sev- 
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also feel the gravitational acceleration from nearby galax- 
ies onto which they are accreted. The question we aim 
to answer is whether or not this additional acceleration 
is the mechanism by which the core of the filament can 
fragment and form dense clouds full of molecules. 

To this end, we have implemented an additional grav- 
itational source term to account for the galaxy. The 
galaxy potential is modeled as an NFW halo (Navarro, 
Frenk, & White 1997) which has a total interior mass as 
a function of radius of 



eral possible mechanisms, including mergers, tidal inter- 
actions, and accretion from surrounding filaments (e.g. 
Cumming et al 2008; James et al 2010). Meier et al 
(2002) also suggest that due to the low metallicity of 
these galaxies, they represent an unique window into the 
early epochs of start formation. 

NGC 5253, a nearby dwarf galaxy at a distance of 3.8 
Mpc (Gibson et al 2000) and companion to M83, is an ex- 
cellent example of such a galaxy. Recent observations of 
this galaxy have uncovered intense bursts of star forma- 
tion. Multiple "super-star clusters" are extremely bright 
in the infrared with luminosities of Ljr ~ 10 9 L© (e.g. 
Gorijian, Turner, & Beck 2001) and visual magnitudes of 
M v ~ -10 to -11 (Meurer et al 1995; Gorjian 1996). Fi- 
nally, dust lanes with motion consistent with inflow have 
been uncovered in radio observations of CO (2-1) (Meier 
et al 2002). However, how much the inflow effects the 
nature of the starburst remains an open question. 

The gravitational potential in the previous simulations 
have all been solely aligned with the radial direction and 
the resulting collapse is fairly uniform. However, in a 
physical situation, this is not the case. These filaments 



M(< R) 



Mtotal 

M totaI 



F(cJ?./fl max ) 
F(c) 



if R < i? m ax, 
if R > -ffmax, 



(24) 



where M to tai is the total mass of the halo, F(t) = ln(l + 
t) — -r^-7 , c is the concentration parameter which we take 
to be 13.5 (Del Popolo & Cardone 2012), and i? ma x = 
kpc is the virial radius of the galaxy 
M to tai/1O 9 M0. This is given as given 



l/3/l+*e 



3^9 v—; 

halo, with Mg 

as 1.14xl0 3 (Lopez-Sanchez et al 2012), where we take 

z c (— 9) to be the virialization redshift. The gravitational 
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acceleration from this profile is then given simply as 



«d 



g,grav 



GM(< R) 

R~ 2 ' 
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Figure [5] shows the ID line profiles for density, H2, and 
CO for both FID and a run with the dwarf galaxy grav- 
itational potential, where profiles in both the r-direction 
(through the center of the simulation, 54.0 kpc away from 
the center of the dwarf halo) and the z-direction (along 
left r-axis). The evolution of the filament remains largely 
the same as the fiducial case. The filament gas recom- 
bines and collapses toward the center of the dark matter 
potential. In both cases the center of the filament be- 
comes very dense and very cold. 

An interesting difference between these two cases, how- 
ever, is the strength of the shock formed as the the fila- 
ment collapses. The strength of the shock is related to 
how quickly the filament collapses. The evolution of the 
fiducial case is outlined in i j4.ll On the other hand, the 
run with the dwarf galaxy the evolution is slightly faster. 
This leads to an increase in the density near the center 
of the filament as compared to the fiducial case, which 
causes the shocks to be less pronounced. However, these 
shocks are largely transient and don't impact the final 
state of the filament in a large way. 

The left panel of Figure [9] shows a 2D logarithmic den- 
sity profile where the center of the galaxy potential is 
at [0,0] at the end of the simulation. The profile shows 
an ideal situation for the continued accretion of gas into 
the dwarf galaxy. The density peaks at the center of the 
dwarf galaxy potential while the rest of the filament gas 
has collapsed into a dense, cold structure stretching along 
the interior r-axis. Once the gas has collapsed into the 



center of the dark matter potential, it is susceptible to 
collapse into the dwarf galaxy. This influx of cold dense 
gas promotes further star formation within the galaxy. 

Finally, the right panel of Figure [9] shows the compar- 
ison of the accreted mass in two equally sized regions 
in the filament. The solid (blue) line shows a [108,108] 
pc region centered on [54,54] pc and includes the re- 
gion nearest the center of the dwarf galaxy. The dashed 
(green) line shows a comparable region centered on [54 
pc, 107.8 kpc] which is likely to be the least altered by 
the gravity from the dwarf galaxy. We plot the evolution 
of the enclosed mass these regions versus time once the 
filament starts to collapse. As expected, the addition of 
the dwarf galaxy drastically increases the amount of gas 
collected in that region. From this, we can estimate an 
accretion rate of ~ 0.1 M Q /yr which, if most of this gas 
is turned to stars, matches the star formation estimate 
given by Lopez-Sanchez et al (2012) (~ 0.2-0.1 M /yr). 

Many observations of NGC 5253 have uncovered multi- 
ple giant molecular clouds (GMC) that are found in and 
around the central starburst and are thought to be the 
product of infalling gas (Meier et al 2002). The forma- 
tion of these GMCs are an open question, although some 
have suggested that they were formed from an interac- 
tion between NGC 5253 and M83 (e.g. van den Berge 
1980; Kobulnicky & Skillman 1995; Calzetti et al 1999). 
Although most of the GMCs are not resolved in their 
radio observations, Meier et al (2002) reports estimated 
physical sizes of these clusters to be < 89 pc in diameter. 
Studies of other, less active, dwarf galaxies have GMCs 
with sizes between ~ 10-70 pc (e.g. Rubio, Lequeux, & 
Boulanger 1993; Wilson 1994; Wilson 1995; Taylor et al 
1999; Walter et al 2001; Meier, Turner, & Beck 2001). 
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To determine whether or not this occurs in our simu- 
lations, we reran the dwarf galaxy simulation with an 
additional two levels of refinement. This gives us a reso- 
lution of 3.5 pc at the finest level. The evolution of the 
filament remains largely the same save the higher central 
density (by about an order of magnitude) at the center 
of filament at the end of the simulation. 

Unfortunately, such clusters are not formed in any 
of the simulations presented. On the other hand, it is 
still unclear to what extent their absence is due to the 
axisymmetric geometry imposed, which may underesti- 
mate fragmentation due to cooling instabilities as the 
filament collapses. Furthermore, due the small sizes of 
some GMCs, even in the extremely-high 3.5 pc resolution 
case, it is possible that we are not resolving the forma- 
tion of these GMCs due in part by the limit on refinement 
level and use of the artificial pressure. Future, extremely 
high-resolution, fully three-dimensional simulations will 
be necessary to settle this issue definitively. 

6. DISCUSSION & CONCLUSION 

Filaments are ubiquitous structures that stretch be- 
tween the galaxies and dark matter halos found in obser- 
vations and large-scale N-body cosmological simulations. 
In galaxies like the Milky Way, they provide the gas reser- 
voir to continue the observed near constant star forma- 
tion. In smaller galaxies, however, the effect is not quite 
as peaceful. Dwarf galaxies experience episodes of star 



formation that are just as strong as those seen in larger 
galaxies. In addition, dwarf galaxies are more metal-poor 
than their larger counterparts, which provides a unique 
window into the process of star formation at earlier star 
formation epochs. However, the mechanism that causes 
such extreme star formation in these small galaxies is not 
well understood, although mergers and infalling gas may 
be viable candidates. 

In this paper we have studied how these initially ion- 
ized, hydrostatic filaments evolve in time by tracking the 
formation and evolution of key molecular coolants and 
their associated cooling. In the fiducial case, much of 
the early evolution is determined by atomic cooling from 
metals and hydrogen-helium lines. As soon as the fila- 
ment gets dense enough important molecules are formed 
which then allows the filament gas to cool to below the 
hydrogen-helium floor. The initial cloud has collapsed 
into the center of the original filament and is character- 
ized by having a very dense and cold core. This trans- 
formation takes nearly a single free-fall time to occur. 

The inclusion of a pressure floor to preserve the Jeans 
criterions prevents us from following the specific details 
of the filament once the gas becomes very cold and dense. 
With improved resolution, the collapse could be contin- 
ued down to much higher densities, to determine in bet- 
ter detail if the conditions and scales are reached at which 
stars are expected to form. However, pushing to such 
high densities is beyond the scope of this paper. 
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profiles. Again, we take the logarithm along the radial direction to highlight the profiles near the center of the filament (Columns 1 and 2) 
or in the center of the dwarf galaxy potential (Column 3) . 



The inclusion of a dissociating UV background has 
little impact on the final state of the filament. As ex- 
pected the inclusion of the background impacts the time 
at which molecular cooling becomes efficient as some 
molecules are destroyed. However, the final state of the 
filament remains the same as the fiducial case. 

This story is not repeated for other initial conditions. 
In both HighT and Highz the filament does not evolve 
drastically from its initial profile. This is due to the 
extremely long cooling time scales encountered during 
the collapse. In the HighT case, the low density and 
high temperatures of the filament results in a cooling 
time scale that can be up to several orders of magnitude 
greater than the free-fall time scale. Highz begins near 
the atomic hydrogen-helium cooling floor, and thus does 
not increase its density sufficiently during the atomic 
cooling stage for molecular coolants to efficiently con- 



tinue the collapse. 

Finally, we studied the effect that the gravitational po- 
tential from a nearby dwarf galaxy has on the collapse of 
the filament. While much of the evolution remains the 
same as the fiducial case, a density peak forms at the 
center of the galaxy potential. Due to the increase in 
density, and hence number densities, this region undergo 
quicker chemical and thermal evolution than the rest of 
the filament. 

We also compared this situation to NGC 5253, a star- 
bursting galaxy located nearby at a distance of 3.8 Mpc. 
The starburst is thought to be powered by infalling fil- 
amentary gas similar to the situation studied, and in 
fact, our simulations reproduce a star formation rate of 
~ 0.1 M /yr if we assume all of the infalling gas is con- 
verted to stars, which is similar to the estimate given 
by Lopez-Sanchez et al (2012). On the other hand, the 
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multiple GMCs seen in radio observations are not repro- 
duced in our simulations. This may be due to the need 
for full three-dimensional simulations and perhaps even 
higher resolutions to fully capture the fragmentation pro- 
cess during the collapse of the filament and its accretion 
onto the dwarf galaxy. For now, how these extragalactic 
GMCs are formed remains an intriguing open question 
to be answered by future studies. 
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